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ABSTRACT 

We consider two kinds of problems: the computation of poly- 
nomial and rational solutions of linear recurrences with coef- 
ficients that are polynomials with integer coefficients; indef- 
inite and definite summation of sequences that are hyperge- 
ometric over the rational numbers. The algorithms for these 
tasks all involve as an intermediate quantity an integer N 
(dispersion or root of an indicial polynomial) that is poten- 
tially exponential in the bit size of their input. Previous 
algorithms have a bit complexity that is at least quadratic 
in N. We revisit them and propose variants that exploit 
the structure of solutions and avoid expanding polynomials 
of degree N. We give two algorithms: a probabilistic one 
that detects the existence or absence of nonzero polynomial 
and rational solutions in 0{\N log 2 N) bit operations; a de- 
terministic one that computes a compact representation of 
the solution in 0(N log 3 N) bit operations. Similar speed- 
ups are obtained in indefinite and definite hypergeometric 
summation. We describe the results of an implementation. 

Categories and Subject Descriptors: 1.1.2 [Symbolic 
and Algebraic Manipulation]: Algorithms 

General Terms: Algorithms, Experimentation, Theory 

Keywords: Computer algebra, polynomial and rational so- 
lutions, linear recurrences, summation, creative telescoping, 
complexity 

1. INTRODUCTION 

A central quantity for many algorithms operating on lin- 
ear recurrences and their solutions is the dispersion. 

Definition 1. The dispersion set of two polynomials P 
and Q in Q[n] is the set of positive integer roots of the re- 
sultant R(h) — Res n (P(n), Q(n + h)). When this set is not 
empty, its maximal element is called the dispersion of P 
and Q. 

Thus, the dispersion is the largest integer difference between 
roots of P and Q. As shown by the simple example (P, Q) = 
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(n, n — N) with N £ N, the dispersion can be exponentially 
large in the bit size of the input polynomials. It cannot get 
much worse: when the polynomials have integer coefficients 
whose absolute value is bounded by B, their dispersion is 
bounded by 4B [12, Fact 7.11]. This exponential size yields 
the dominant term in the worst-case complexity of many 
algorithms computing — or operating on — solutions of 
linear recurrences. 

For instance, the computation of a Gosper-Petkovsek form 
produces a polynomial whose degree N can be linear in the 
dispersion of its input and thus exponential in its bit size. 
If this polynomial is expanded it has N + 1 coefficients; over 
the integers, its total bit size is O (N 2 log N). This form is 
used in the first step of Gosper's summation algorithm and 
of Abramov's algorithm for computing rational solutions of 
linear recurrences. Thus, it makes an important contribu- 
tion to the complexity of these algorithms. Once this form is 
computed, these algorithms search for polynomial solutions 
of an associated linear recurrence. This is done by linear 
algebra using a bound on the possible degree of solutions 
which is at least as large as N, leading again to a more than 
quadratic complexity, even when no nonzero solution ex- 
ists. In turn, a parameterized variant of Gosper's algorithm 
forms the basis of Zeilberger's definite summation algorithm 
which inherits this costly behaviour. By contrast, we pro- 
vide a probabilistic algorithm that detects that no nonzero 
rational solution of a homogeneous linear recurrence exists 
in (D(y/~N log 2 N) bit operations and a deterministic algo- 
rithm that gives a compact representation of all solutions 
in 0(N log 3 N) bit operations. All the algorithms in the 
present work eventually rely on the computation of polyno- 
mial solutions of linear recurrences. In a previous work [7], 
we dealt with the analogous problem in the linear differ- 
ential case, by exploiting the linear recurrence satisfied by 
the coefficients of power series solutions and reducing the 
computation to that of matrix factorials. For the latter op- 
eration, there exist fast probabilistic and deterministic al- 
gorithms (see [8, 10] and the references in [7]). In the case 
of linear recurrences, it is not true that the coefficients of 
polynomial solutions satisfy a linear recurrence in general; 
however, it becomes true if the polynomials are expanded 
in a binomial basis [5, Ch. XIII, art. 5]. This is the basis 
of a simple quadratic algorithm [3] to compute polynomial 
solutions. In Section 2, we show how this conversion is per- 
formed, we recall the basic results on matrix factorials and 
apply them to get the announced complexities. 

From there, in Section 3, we proceed in three steps: (i) we 
slightly modify the computation of the Gosper-Petkovsek 



form so that it does not expand the potentially large poly- 
nomial but instead computes a first-order, moderately-sized 
recurrence for it; (ii) we show that this first-order recurrence 
can be used to compute a linear recurrence satisfied by the 
numerators of rational solutions, in a complexity that is only 
logarithmic in N, both in the homogeneous and nonhomoge- 
neous cases; (Hi) we then compute the numerators as poly- 
nomial solutions via matrix factorials. The close relation 
between Abramov's and Gosper's algorithm makes it possi- 
ble to transfer these results to Gosper's algorithm. Then in 
Section 4, we show how this machinery can be adapted to 
the parameterized variant needed in Zeilberger's algorithm. 
Finally, we describe experimental results in Section 5. 

Notations and complexity measures. All along this 
text, 7Z denotes a linear difference operator with coefficients 
in Z[n]. We view it as a polynomial in the non-commutative 
ring Q(n, S n ), where S„ is the shift operator S n u(n) = u(n+ 
1). Similarly, S x , Sk, and S m denote the shifts with respect 
to x, k, and m. To any difference operator 7Z is attached a 
homogeneous linear recurrence equation TZu = 0. We view 
the solution u either as a sequence (u n ) (also denoted u n ), 
or as a function u(n) (the cases of particular interest being 
polynomial and rational functions). 

For our complexity analyses, the measure we use is the bit 
(or boolean) complexity. For this purpose, our complexity 
model is the multi-tape Turing machine, see for instance [f 8] . 
We use the number of bit operations to express time com- 
plexities in this model. We call bit size (or simply size) of an 
integer a/0 the number X(a) := [log |a|J + f (logi denotes 
the logarithm of x in base 2). By convention, we assume 
that A(0) = f . The bit size of a matrix or vector is the sum 
of the bit sizes of its entries. Polynomials given as input to 
our algorithms are stored in a dense representation; a mea- 
sure of their bit size is given by the sum of the bit sizes of 
their coefficients, including the zero coefficients. Similarly, 
the bit size of a linear recurrence equation (LRE) is the sum 
of the bit sizes of its polynomial coefficients. 

To simplify complexity estimates, we assume that the 
product of two integers of bit size d can be computed within 
1(d) = C(dlog dlog log d) bit operations using Fast Fourier 
Transform [19]. To keep the notation compact, we some- 
times write 1(d) = (5(d); the tilde indicates that the factors 
polynomial in log d or smaller have been omitted. 

For any prime number p, the bit complexities of the op- 
erations (+,— , X,-r) in the finite field F p := Z/pZ are in 
0(1 (log p) log log p). We assume that over the rings we use, 
the product of two polynomials of degree at most d can be 
computed within 0(M(d)) base ring operations (each ring 
operation being counted at unit cost) and that M(d) = 
0(d) [9, ch. 2]. For computations in F p [x], the bit complex- 
ity is bounded by multiplying the arithmetic cost estimates 
by the bit complexity of the basic operations in F p . 

In all our algorithms we are interested in reducing the 
complexity with respect to a potentially exponential param- 
eter x (related to a dispersion or to a root of an indicial 
polynomial). Thus we consider as having cost 0(1) any op- 
eration whose complexity is polynomial in the bit size of the 
input recurrence or polynomials, and concentrate on the de- 
pendency of the complexity in i. In order to provide the 
code with an actual bound on the size of primes that need 
to be used so that the bound on probability of error is guar- 
anteed, we have to perform a much more precise complexity 



analysis taking into account all parameters (order, degree of 
coefficients) (as in the proof of [7, Thm. 3]). Such a detailed 
analysis will appear in [6]. 

2. POLYNOMIAL SOLUTIONS 

In symbolic summation and in the resolution of linear re- 
currences, all the known algorithms ultimately require poly- 
nomial solutions of linear recurrence equations. 

In this section, we give algorithms for computing descrip- 
tions of the Q- vector space of solutions of a linear recurrence 
operator 1Z with coefficients in Z[n]: 

TZu = a r (n)u(n + r) + • • • + ao(n)u(n) =0, n > 0. (1) 

We focus on two types of solutions of such recurrences: so- 
lutions with finite support and polynomial solutions. 

In what follows, we make the hypothesis that is an or- 
dinary point of the recurrence. This means that the lead- 
ing coefficient a r (n) does not vanish at any of the integers 
0, 1,2,.. .; in other words, when unwinding the recurrence, 
no division by zero is encountered. This condition is ensured 
after a generic translation hh n + a. Under our complexity 
assumptions, a proper a and the corresponding translation 
can be computed in a polynomial number of bit operations, 
so that there is no loss of generality for the problems we 
consider. The general case (when is not ordinary) is tech- 
nically more demanding but does not change the complexity 
estimates we give here. It will be presented in [6]. 

Let Sol(7?.) denote the vector space of solutions (uo, Ui, ■ ■ .) 
of (1). In the case of an arbitrary 7Z, the dimension of 
Sol(72.) as a Q-vector space may be different from r (both 
larger or smaller). However, when is an ordinary point 
of 7Z, Sol(72.) has dimension exactly r and a basis is given by 
the sequences it"' , j = 0, . . . , r — 1, satisfying 1Z and having 
initial conditions — for < i < r— 1, where S m ,n 

is Kronecker's 5 symbol (<5 m , m = 1, S m , n — if m 7^ n). 

In §2.1, we describe the compact representation that forms 
the basic data structure of our algorithms. Then, in §2.2 
we recall classical results that allow for the efficient com- 
putation of the Nth element of a solution of 1Z. In §2.3 
we describe the reduction from the problem of searching for 
polynomial solutions to that of finding solutions with finite 
support. Next, we give in §2.4 algorithms to compute finitely 
supported and polynomial solutions of recurrences. We con- 
clude this section by showing in §2.5 how the evaluation of 
a polynomial and its finite differences can be performed ef- 
ficiently in the compact representation. 

2.1 Compact Representation 

Classically, a polynomial solution u(n) of (1) is repre- 
sented by its coefficients in the monomial basis {n k }. We use 
an alternative data structure for u(n), which is motivated 
by the observation that its coefficients ct in the binomial 
basis {(£)} obey a recurrence with polynomial coefficients. 

Example 1. The recurrence (n + l)u(n + 1) — (n + TV + 
l)u(n) = has a unique nontrivial monic polynomial solu- 
tion u(n) = (n+1) • • • (n+N). To write down its coefficients 
in the monomial basis at least ^N 2 logN bits are needed. In 
contrast, u{n) can be represented by the recurrence 

(k + l)cfc+i - (N - k)c k = 0, c =N\ 

on the coefficients Ck ofu(n) in the binomial basis {(£)}; the 



bit size of this new representation is only linear in NlogN, 
most of the size being in the initial condition Co = AH. 

Definition 2. The compact representation of a polyno- 
mial solution of (1) is the data of a linear recurrence and 
initial conditions for its coefficients in the binomial basis, 
together with an upper bound on its degree. 

Our aim in this article is to demonstrate that this represen- 
tation of polynomial solutions of recurrences can be carried 
through different algorithms from indefinite and definite hy- 
pergeometric summation and that it is beneficial from the 
complexity point of view. The reason why this representa- 
tion deserves the name "compact" appears in §2.4 below. 

2.2 High-Order Terms of Sequences 

Let (w„) be a sequence satisfying (1). The recurrence TZ 
can be rewritten as a first-order matrix recurrence U n +i = 
C(n + l)U n , where U„ is the vector {u n ,u n -i, ■ ■ ■ , «»-r+i) 
and C is an r x r matrix with rational function entries. 
The problem of computing a selected term um reduces to 
that of computing U r -i and the matrix factorial F(N) := 
C(N) ■ ■ - C(r). This makes sense since under our hypothesis 
the leading term of the initial recurrence does not vanish 
at 1,2, ... ,N. The numerator and denominator of the ma- 
trix factorial can be computed efficiently, either in Z using a 
binary splitting algorithm, or modulo a prime p using a baby- 
step/giant-step algorithm. These algorithms are described 
in [7, §2.1, §3.1], see the references therein. For further use, 
we extract from [7] the following result. 

Theorem 1 ([7]). Let (Ui) be a sequence of vectors of 
rational numbers that satisfies a recurrence Ui+i — C(i + 
l)Ui, withC(x) anrxr matrix with rational function entries 
in Q(x). Letp be a prime number such that the denominator 
of C does not vanish modp at 1, 2, . . . , N. Then, as N — > oo: 

(a) T{N) = C(JV) • • • C(r) andU N have bit size 0(N log N); 
their values can be computed using O (I (AT log N) log AM 
bit operations. 

(b) J-(N) modp and Un modp can be computed using 

bit operations. 

(c) The rank of J-(N) can be computed in (D(\/~N) bit op- 
erations using a probabilistic Monte Carlo algorithm. 

2.3 Expansion in the Binomial Basis 

For completeness, we recall here an algorithm from [4] 
that we call RecToRec to perform the conversion from a re- 
currence with polynomial coefficients to the recurrence sat- 
isfied by the coefficients of series solutions in the binomial 
basis. Earlier (and slightly more complicated) algorithms 
have been given in [5, Chapter XIII] and [3, Section 4.2]. 
The starting point are the following two identities: 
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If u(n) = J21 



|Cfc 



then applying these identities to 



rewrite u(n+l) and nu(n) and extracting coefficients of (?) 
shows that the ring morphism <j> : Q[n, S n ] — > Q[k, Sk, S^ 1 } 
defined by 4>(S n ) — 1 + Sk and 4>{n) = k(l + S^ 1 ) sends a 



homogeneous LRE satisfied by u(n) to another one satisfied 
by Cfe. The image of (1) is a LRE of the form 

(a r (k)S r k + br-iifyS^- 1 + ■■■ + b- s (k)S^ 3 )c k = 0, k > 0, 

where the leading term is exactly that of (1) and the trailing 
term may involve a negative shift (by convention, c k — 
when k < 0). In particular, if is an ordinary point for TZ, 
so is it for (f>(TZ). The resulting algorithm is as follows. Its 
complexity is clearly polynomial in the bit size of TZ. 



Algorithm RecToRec 

Input: a recurrence TZu = 0, where u(n) — J2k ^{Z)- 
Output: a recurrence S satisfied by the sequence (ck), 
plus a set £ of linear equations on its initial conditions. 

1. Compute T = 4>(TZ); 

2. Let —s — vals fc (T) be its valuation w.r.t. Sk', 

3. If s < return S :=T and £ := 0, 

4. Otherwise return 5 := S^T and the equations £ := 
{((SiT)c)\ k=o = 0,i = 0,...,s-l}. 



2.4 Finite Support and Polynomial Solutions 

We consider here the problem of computing a basis of 
solutions with finite support, that is, whose terms beyond a 
certain index are all zero. The degree of a solution with finite 
support u is, by definition, the unique integer n such that 
u n 7^ and u n +» = 0, for all i > 1. A universal bound N on 
the degrees of all solutions with finite support of the input 
recurrence TZu = is given by the largest positive integer 
root of the trailing coefficient ao(n) of 1Z. Note that N is 
generally not bounded polynomially in the bit size of 7Z. 

Recall that Sol(7?.) has dimension r, with a basis B formed 
by the sequences u^\ with < j < r — 1, given by the ini- 
tial conditions up = 5j,i, for < i < r — 1. Thus, a 
finitely supported solution u is (an unknown) linear com- 
bination 5ZJ=o ^3 U ^^ such that the elements in the slice 
un+i, ■ ■ ■ , UN+r all vanish. This yields linear constraints 
on the initial conditions Aj . 

To determine these constraints, it is sufficient to compute 
the values at indices N + 1 , . . . , AT + r of all the elements 
in B using Thm. 1. The rank of the resulting r x r matrix 
gives the dimension of the vector space of solutions with 
finite support. Since the entries of this matrix have bit size 
0(N log N), the desired Aj's, which are determined by a 
kernel computation, also have bit size 0(N log N). Putting 
together these considerations, we get the following result. 



Theorem 2. There exists a basis (u 



id 



Ad) 



) (d < 



r) of solutions of recurrence (1) with finite support, where 
each Vj is uniquely specified by the data of initial condi- 
tions of bit size 0(N log N) , with N a bound on the inte- 
ger roots of ao . The dimension d as well as the degrees of 
the 's can be computed by a probabilistic algorithm using 
0(M(v / A : )l(log N)) bit operations. The initial conditions of 
the uS 1 ' 's and their maximal degree D < N can be computed 
deterministically in 0(\(Dlog D) logD) bit operations. 

Thm. 2 is the basis for using the name "compact representa- 
tion" : it shows that the compact representation has a size of 
the same order as the initial conditions, while the expanded 



polynomials have size 0(N 2 log N). In general, this latter 
bound is reached. 

Using the results in §2.3, Thm. 2 carries over literally to 
the compact representation of a basis of polynomial solu- 
tions of the recurrence IZu = 0. The corresponding state- 
ment requires a bound on the degree of polynomial solutions 
that is given by the roots of the indicial polynomial. 

Definition 3. The indicial polynomial of 1Z at infinity 
is the trailing coefficient of RecToRec(72.) . 

Corollary 1. The statement of Thm. 2 holds for poly- 
nomial solutions of IZu = 0, with N the largest integer root 
of the indicial polynomial of 1Z at infinity. 

Nonhomogeneous Equations. We now consider the equa- 
tion lZu(n) — f(n), with coefficients in Z[n] and right-hand 
side of degree m. Applying RecToRec and expanding f(n) 
in the binomial basis, the initial problem boils down to 
the search of finitely supported solutions of a nonhomoge- 
neous equation Sc(k) — g(k), where g is a sequence with 
finite support, g(i) — for i > m. In matrix notation, 
we have Uk+i = C(k + l)Uk + ffc+i, where Uk is the vector 
(tifc, . . . , Uk-r+iY and v k is the vector (g(k),0, . . . , 0)'. Then 
the vector of initial conditions U r -i satisfies the affine con- 
straint A(BU r -i + w m ) = 0, where A := C(N + r)C(N + 
r - 1) • • • C{m + 1), B := C(m)C(m - 1) • • • C(r), w, := 
Vi + C(i)wi-i, for r + l<i<m + l and w r — v r . 

For large N, using Thm. 1, the matrices A and B can be 
computed efficiently. The bit size and the computational 
cost of ™ m +i is 0(1). Thus, solving the affine system of 
size O(N) yields the finitely supported solutions of Sc — g 
and the polynomial solutions of IZu = / and we get the 
following. 

Corollary 2. Let f be a polynomial. Then the state- 
ment of Thm. 2 holds for polynomial solutions of nonho- 
mogeneous equations TZu(n) — f(n) as the largest integer 
root N of the indicial polynomial of 1Z at infinity becomes 
large. 

2.5 Evaluation in Compact Representation 

The compact representation is not only a data structure 
for intermediate computations. It can actually be exploited 
further. In particular, we now detail the evaluation at an 
algebraic number a of a polynomial u(x) and an iterated 
difference A H (u) (where A = S x — 1 and H is potentially 
large). The polynomial u is given by its degree iV and the 
recurrence 

r 

^2 ai(k)c(k + i) =0 forallfc>0 

satisfied by its coefficients Ck in the binomial basis {(?)}, 
together with initial conditions. The basic idea is embodied 
in the following. 

Lemma 1 (Folklore). If (uk) and (vk) are solutions 
of linear difference equations with polynomial coefficients, 
then so is the sequence (ujv) defined by un = ~^2,k=o UkVk - 

This lemma can be applied to the sequences (c&) and (^). 
Evaluating the resulting sequence at N using Thm. I gives 
u(a) for 0(\(N log N) log N) bit operations, when N is large. 



Using Pascal's formula A H = f . , we deduce that 

A H u(a) = X^fc^Tt/ 7 Ck + H i'k)' recurrence satisfied by the 
sequence {ck+n)k is obtained by shifting by H the coeffi- 
cients of the recurrence of (ck). This new recurrence has bit 
size 0(logH) and initial conditions can be determined by 
binary splitting in 0(1 (N log N log H) log N log H) bit oper- 
ations. Here, our asymptotic bound involves the two pa- 
rameters N and H, as both are potentially exponential in 
the input size. As above, the compact representation of the 
recurrence satisfied by Dk '■= J2e=o °i+h (") can be deter- 
mined efficiently, as well as its iVth term A H u(a). 

3. RATIONAL SOLUTIONS 

3.1 Compact Gosper-Petkovsek Normal Form 

The classical Gosper-Petkovsek normal form [16, 14] of 
a reduced rational function P/Q in Q(n) consists of three 
polynomials A, B, C in Q[n] such that 

P(n) _ A(n) C(n+1) 

Q(n) B(n) C(n) ' 1 ' 

with the constraints 

gcd(A(n), C(n)) = 1, gcd(B(n), C(n + 1)) = 1, 

and for all h € N, gcd(A(n),B(n + h)) = 1. (3) 

The degree N of the polynomial C(n) is potentially expo- 
nentially large. Thus, in our algorithm CompactGPF below, 
we modify the usual algorithm (e.g., in [17]) slightly so that 
the polynomial C(n) is not expanded. Similar ideas appear 
in [13] in the context of indefinite rational summation. 



Example 2. CompactGPF(n, n - N) = (1, 1, {(n, N)}). 

Note that the input is an ordered pair (P, Q) and not a ra- 
tional function P/Q. The output of the algorithm changes 
if (P, Q) is replaced by (FP, FQ) for F € Q[n]. This will be 
necessary for our treatment of rational solutions below. On 
the other hand, the output A, B, and (/;'s also satisfy (3) 
whenever P and Q have no common factor, so that the 
Gosper-Petkovsek normal form of a rational function in Q(n) 
given in reduced form P/Q is obtained by CompactGPF(P, Q) 
As an outcome of this algorithm, the rational function 
C(n) / C(n + j) (j = 1, 2, . . . ) is easily obtained as 

C{n) = A gi{n + j-l-hi)---gi(n-hj) 
C(n + j) f = l gi ( n + j-l)... 9i (n) ■ {> 



Algorithm CompactGPF 

Input: an ordered pair (P(n),Q(n)) of polynomials. 
Output: (A(n),B(n),{(gi(n),hi),i — l,...,s}) such 
that C(n) — Yli gi( n — 1) ' ' ' gi( n — hi) satisfies (2). 

1. Compute hi > • • • > h a > the positive integer 
roots of ReSn(P(fi), Q(n + h)); 

2. A(n) := P{n),B{n) := Q(n); 

3. For i from 1 to s do 

a. gi(n) := gcd(A(n),B(n + hi)); 

b. A(n) := A(n)/ gi (n), B(n) := B(n)/ gi (n - hi); 

4. Return (A, B, {(ffi(n), hi),i = 1, . . . , s}). 



Algorithm HomCompactRatSols 

Input: a homogeneous LRE lZu(n) — 0. 

Output: a basis of its rational solutions in compact form 

1. (A, B, C) := CompactGPF(a r (n - r + 1), a (n)); 

2. Normalize C(n)lZ(v(n)/C(n)) using (4) and denote 
the result Tv(n); 

3. Compute a basis B of the polynomial solutions of 
Tv(n) = 0; 

4. Return {p(n)/C(n) \ p(n) £ B}. 



For large iV and j — 0(1), it has "small" numerator and de- 
nominator of degrees bounded by j times those of P and Q. 
This equation for j = 1 is a homogeneous LRE that plays the 
role of a compact representation of C. The initial value C(0) 
(more generally C(k) where k = 0(1)) has size 0(N log N) 
and can be computed by Thm. 1 within 0(I(N log N) log N) 
bit operations. In the next sections, we use this to design 
"compact" variants of Abramov's and Gosper's algorithms. 

Proposition 1. Algorithm CompactGPF is correct. For 
(P, Q) with rational coefficients, it has deterministic polyno- 
mial bit complexity in the bit size of (P, Q) . 

Proof. The correctness is that of the classical algorithm 
since the only difference is that we do not expand C. Step 1 
is dealt with by a deterministic algorithm due to Loos [15] 
(cf. [12, 13] for faster probabilistic algorithms). Step 3 is 
performed at most degPdegQ times, and each step is poly- 
nomial by the classical algorithms as found in [11]. □ 

3.2 Compact Rational Solutions 

We now consider rational solutions of the LRE lZu(n) — 
f(n), with / a polynomial in Q[n]. 

Our starting point is the following result of Abramov [1]. 

Lemma 2 (Abramov). The polynomial C(n) of the 
Gosper-Petkovsek form of (a r (n — r + l),ao(ra)) is a multi- 
ple of the denominator of all rational solutions of lZu{n) — 
fin). 

Abramov's algorithm first computes C(n), then performs 
the change of variable u{n) = v(n)/C(n), leading to 



v(n + r) v(n) 
ar ( n ) 7v — I — T H h ao (n) -— ■ = f{n) , 



C{n + r) 



C{n) 



(5) 



whose polynomial solutions v(n) are then sought. 

In the homogeneous case (f{n) = 0), using (4) reduces 
this equation to an equation of polynomial size. This is 
described in Algorithm HomCompactRatSols (see Figure). 
In Step 2, the "Normalize" operation consists in expand- 
ing C(n)/C(n + j) using (4) and taking the numerator of 
the resulting expression. Also, if necessary, we change n 
into n + a with C(a) ^ 0, so that is not a singular point 
in Step 3. This can be detected and changed at a cost 
of 0(\(N\og N) log N) operations. In Step 4, the output 
is given by the compact forms of the numerators and C is 
given by the output of CompactGPF. 

In the nonhomogeneous case, reducing (5) to the same 
denominator would lead to an equation whose right-hand 



side has a potentially exponential degree. Instead, we con- 
sider the homogeneous operator 5 = (f(n)S n — f(n + 1))TZ, 
whose bit size is polynomial in that of lZu(n) — f(n) and 
that can be treated by the algorithm above. If u n is a ra- 
tional solution of S, then w n = lZu n is a rational solution 
of f(n)w„+i — f(n + l)w n . This implies that w n = A/(n) 
for all n larger than the largest root of / and since w n is 
rational, also for all other values of n. Thus fixing A so that 
TZu(k) — f(k) for any k such that f(k) ^ concludes the 
computation. This is the basis of the following algorithm. 



Algorithm NonhomCompactRatSols 

Input: a LRE TZu{n) = f{n), with / / 0. 
Output: a particular rational solution p and a ba- 
sis (bi,...,bd) of rational solutions of IZu in compact 
form 

1. W := HomCompactRatSols((/(n)5 n - f(n + 1))TZ); 

2. Find k G N such that f(k) / 0; 

3. Write Tl{^2, w&w iww{k)) =: U{£) for an unknown 
£; = {£,w)wew and solve U{£) = for a basis 
(fJr 1 , ■ ■ ■ , A* ) of its solution space and W(£) = 
f(k) for a particular solution A; 

4. Return p := ~^2 wew \ w w{n) and the bi's given by 
6 * ~ Y, w ew a ™ w ( n )- 



In Step 2, just iterating k — 0, 1, . . . till a point where / 
is found to be nonzero is sufficient for our purpose. If 
is a bound on the degree of the numerators and denomi- 
nator computed in Step 1, then the values of the w(k)'s in 
Step 3 have size 0(N log AQ and can be computed by binary 
splitting. From there, it follows that the affine equation in 
Step 3 has coefficients of size 0(N log A^), which is then also 
a bound on the size of its solutions. These solutions can be 
computed in the form of a point and a basis of a vector space 
within 0(\(N log N) log A^) bit operations by standard linear 
algebra. The same complexity is sufficient for the products 
of initial conditions in Step 4. 

The results of this section are summarized as follows. 

Theorem 3. Let N be the sum of the largest nonnega- 
tive integer root of the indicial polynomial of 1Z at infinity 
and the degree of the polynomial C(n) of (2) with P{n) = 
a r (n — r + l) and Q(n) = ao(n). The dimension of the affine 
space of rational solutions oflZu{n) — f{n) can be computed 
probabilistically using O (M(\/~N)\(\og N)) bit operations. A 
compact representation of the solutions can be computed de- 
terministically in 0(\(N log A^) log N) bit operations. 

Proof. The largest integer root of the indicial polyno- 
mial of TL at infinity is a bound on the valuations of power 
series solutions of IZu — at infinity, including the valuation 
of v(n)/C(n). Adding the degree of C gives the announced 
bound on the degree of polynomial w's. From there, the 
theorem follows from Cor. 1. □ 

3.3 A Compact Gosper Algorithm 

Given a hypergeometric term t(n), i.e., such that t(n + 
l)/t(n) =: r(n) £ Q(n), Gosper's algorithm [14] finds its 
indefinite hypergeometric sum, if it exists. Such a sum 
is necessarily of the form u(n)t(n) for some u(n) G Q(n). 



Zeilberger's Algorithm 

Input: two functions f ("+ 1 ' r ") anc j t(n,m+i) j Q(n,m). 
Output: a LRE £[ =0 ^.(™)4(E„ t(n,m)) = 0. 
For r = 0, 1, 2, . . . do 

1. Construct the equation (-E r ) 

t(n + 1, m) t(n,m + i) 

u(n+l,m)—- u{n,m) = } Xi{m)—- — , 

t[n,m) r - ?, t[n,m) 

2. Find if there exist Ai's in Q(m) so that (-Er) admits 
a solution u(n,m) £ Q(n,m); 

3. If so, compute and return them; otherwise proceed 
to the next r. 



Small Linear System 

Input: the equation (E r ) from Zeilberger's algorithm. 
Output: an equivalent system linear in the Ai. 

1. Compute lZu{n) = f(n), the numerator of (E r ); 

2. Compute a multiple C(n) of the denominator of its 
rational solutions and a bound N on the degree in n 
of their numerators; 

3. Compute Sv(n), the numerator of 
C(n)(f(n)S n - f(n + l))1Z{v / C){n); 

4. Compute (T,£) := RecToRec(S); set £ := £ U 
{TZ(v/C)(0) = /(0)}; let a be the order of T; 

5. Compute the value (cn+i, ■ ■ ■ ,cn+s) =: V for a 
nonzero sequence solution of RecToRec(C7?.(u/C)); 

6. Compute the value W := (djv+i, . . . , djv+ s ) for an 
arbitrary sequence solution of T obeying £ ; W is 
of the form W* + Yn=o OIU y ^* depends on 
the initial conditions; 

7. The system (E) := (m^ + ELo = °) is simul " 
taneously linear in the Ai's and /j,. 



Thus, the problem is reduced to finding rational solutions 
of u(n + l)r(n) — u(n) = 1. This can be solved by Non- 
homCompactRatSols. A further optimization is present in 
Gosper's algorithm: if r(n) = P(n)/Q(n) in reduced form, 
the polynomial B(n) of (2) satisfies (3), so that it divides 
the numerator of u(n + 1). (This can be generalized to de- 
tect factors of numerators in arbitrary LRE's). This does 
not affect the expression of the complexity result, which is 
as follows. 

Theorem 4. Let t(n) be a hypergeometric term such that 
t(n+l)/t(n) =: P(n)/Q(n) € Q(n), with gcd(P, Q) = 1. Let 
N be a bound on the degree ofC in (2) and on the largest pos- 
itive integer root of the indicial polynomial of P(n) S„ — Q(n) 
at infinity. Then the existence of an indefinite hypergeomet- 
ric sum of t(n) can be determined by a probabilistic algo- 
rithm using 0(M(\ /r N)I(log N)) bit operations, a compact 
representation of it can be computed deterministically us- 
ing £)(l(JVlog N) log N) bit operations. 

Note that in the special case of rational summation (i.e., 
t(n) £ Q(n)), it is actually possible to decide the existence 
of a rational sum in only polynomial complexity, see [13]. 

4. DEFINITE HYPERGEOMETRIC SUMS 

A bivariate hypergeometric term t(n, m) is such that both 
t(n+l, m)/t(n, m) and t(n, m+1) /t(n, m) belong to Q(n, m). 
Given such a term, Zeilberger's algorithm [21] computes a 
LRE satisfied by T(m) — ^ n t(n, m). The idea is to synthe- 
size a telescoping recurrence, i.e., a rational function u(n, m) 
and a linear operator P(m, 5 m ) such that 

(S n — l)u(n,m)t(n,m) — P(m, S m )t{n,m). 

Indeed, summing over n and granted boundary conditions 
known as "natural boundaries" , we obtain P(m, Sm)T(m) = 
0. If P was known, then Gosper's algorithm would find the 
left-hand side. This is the basis of Zeilberger's algorithm 
(see Figure). Termination is guaranteed only if such a LRE 
exists. This occurs in the so-called "proper-hypergeometric" 
case [20] and a general criterion has been given by Abramov [2] 

Note that knowing it permits to check the output oper- 
ator P by simple rational function manipulations, which is 
why the rational function u is called "certificate" in [17]. 

Zeilberger's algorithm is based on a refinement of Gosper's 
algorithm for Steps 2 and 3. It reduces the computation in 



Step 2 to solving a system that is linear simultaneously in 
the Ai's and in another set of A'"-!- 1 variables, where A^ is po- 
tentially exponential in the bit size of (E r ), see e.g. [17, §6.3]. 
An equivalent linear system in a small number of variables 
can be computed by Algorithm Small Linear System (see Fig- 
ure) . The important point is linearity: not all solutions of T 
are linear in the Ai's, but this property is ensured when the 
initial conditions satisfy £. Indeed, in Step 2, by Lemma 2, 
C does not depend on the Ai's. Then, by induction on n, 
starting from lZ(v/C)(0) — /(0), the factor f(n) of the lead- 
ing coefficient in 5 cancels out and thus the solution v(n) 
is linear in the A,:. This property is then preserved by the 
linearity of RecToRec. The final system (E) has solutions if 
and only if (E r ) has rational solutions. 

The description of Small Linear System is geared towards 
the use of compact representations and matrix factorials in 
intermediate steps. This is straightforward for Steps 1-5. In 
Step 6, we cannot make direct use of the factorial of the ma- 
trix associated to T: this matrix involves the Ai's rationally 
and its factorial has too large a size for our target complex- 
ity. Instead, we exploit the linearity in the Ai's by construct- 
ing the vector W using matrix factorials for A a vector of 
0's with a 1 in zth position for i — 0, . . . ,r and setting the 
initial condition to 0, which gives the coefficients Wi. 

From there we derive our compact version of Zeilberger's 
algorithm given in Compact Zeilberger Algorithm. In Step 2, 
the whole construction can be performed by matrix factori- 
als with integer entries, within the complexities of Thm. 1. 
If a rational solution (Ai(m)) exists, then the system (E) 
has the corresponding (Ai(mo)) for solutions. Thus if (E) 
does not have a nonzero solution, (E r ) does not have a ra- 
tional one. This gives a fast probabilistic test. Then, in 
Step 5, the algorithm is used again with matrices that are 
polynomial in the variable m. In that case, the system (E) 
can be computed by binary splitting with 0(M(N) log A^) 



Compact Zeilberger Algorithm 

Input: two functions Us+hlsH anc j t(n,m+i) ■ Q(n,m). 
Output: a LRE £Lo AiMS^Q^ *0, m)) = 0. 
For r = 0, 1, 2, . . . do 

1. Take a random mo £ Q and construct (E r ) 
with m = mo; 

2. Apply Small Linear System to this equation; 

3. Find if there exist nonzero solutions to this system; 

4. If not, proceed to the next r; 

5. Otherwise, construct (E r ), apply Small Linear Sys- 
tem and return its solutions. If it does not have 
nonzero rational solutions, go to Step 1. 



arithmetic operations. The final system has coefficients of 
degree 0{N) with coefficients of bit size <D(N log N) each 
and this is also the size of the Ai's to be found. At the same 
time, we find fi, which gives us a compact representation of 
the certificate. 

An optimization is obtained by using the values of the 
Ai(mo)'s to compute the value N' of the degree of the corre- 
sponding sequence. With high probability this is the actual 
degree in n of the numerator of u(n, m), which can be much 
smaller than N, thus saving a lot of computation in Step 5. 

The following theorem summarizes this section. 

Theorem 5. Let t(n,m) be hypergeometric over Q. Let N 
be the maximal number of variables in the linear system 
solved in the classical version of Zeilberger 's algorithm. Then 
it is possible to detect probabilistically that this system does 
not have any nonzero solution in 0(M(\/ r N)\(log N)) bit op- 
erations. If it does have a solution, it is possible to com- 
pute the corresponding Xi 's of degree O(N) and total bit size 
0(N 2 logN), as well as a compact representation of the cer- 
tificate, in 0(M{N)\(N log N)) bit operations. 

For the sake of comparison, a crude analysis by unrolling the 
triangular system of dimension N + r + 2 and taking into 
account coefficient growth leads to a 0(N 4 ) bit complexity 
estimate for the classical algorithm, which can be reduced 
to 6(N S ) by using the binomial basis. 

5. EXPERIMENTAL RESULTS 
5.1 Rational Solutions 

We consider two families of linear recurrences: 

2n{N - n)(-4N - 3nN + 6 + 3n 2 + 8n)u(n) 
- (n + l)(-3nN + 2N + 3n 2 - An- 4)(n + 1 - N)u(n + 1) 
+ (n + 2)(-3n7V - N + 3n 2 + 2n + l)(n + 2- N)u(n + 2) = 0, 

2n{n - 2N)(n - N)(n 2 - 3nN + 3n + 27V 2 - 3N + 2)u(n) 
-(n+l)(n+l-2N)(n+l-N)(3n 2 +6n-9nN+6N 2 -4N)u{n+l) 
+ {n+2)(n+2-2N)(n+2-N){n 2 +n-3nN+2N 2 )u(n+2) = 0. 

The first one (-Ri) does not have any rational solution, while 
the second one (ife) has l/(n(n — 2N)) as a solution. In both 
cases, when N is a large integer, a large dispersion has to 
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Table 1: Timings (in sec.) for classical and compact 
versions of Abramov's algorithm 

be considered. In Table 1, we give a comparison of the tim- 
ings obtained by our Maple prototype (denoted Compact) 
and that of the command ratpolysols of Maple's pack- 
age LREtools (denoted Classical). This table illustrates the 
"nonexponential" character of the compact versions of the 
algorithms. In the first case, both output are identical (no 
solution). In the second case, however, we return a compact 
representation of the output. For instance, with N = 2 100 
we get (in 0.04s) the denominator (n(n- 2 100 )(n - 2 101 )) (in 
expanded form) and for the numerator the recurrence 

(1 - k 2 )c k + (2 100 + Ak - k 2 )c k+1 

+ {k 2 -2k- B)c k+2 + (k 2 -Ck + D)c h+3 = 0, 

satisfied by its coefficients in the binomial basis, together 
with initial conditions Co = — 2 100 , c\ = 1, where the coeffi- 
cients A, . . . ,D are 200 bit long integers. 

5.2 Definite Hypergeometric Summation 

We consider the following family of hypergeometric terms: 

(2n + m + N\ ( 2m\ / m\ 
^™)=[ N )\2n){n)- 

For /YeN, the sum ^2 n t(n, m) satisfies a third-order homo- 
geneous LRE. When Zeilberger's algorithm is executed on 
this term, the bound it has to use on the degrees of numera- 
tors of rational solutions of the equation (E r ) is JV + 3(r — 1). 
This plays the role of a "large" N and makes it possible to 
exhibit the complexity behaviour of the algorithms. 

In Table 2, we give a comparison of the timings obtained 
by our prototype implementation in Maple (denoted "Com- 
pact") and those obtained by Maple's Zeilberger command 
in the package SumTools : -Hypergeometric (denoted "Clas- 
sical"). The indication "> 2Gb" means that the compu- 
tation had to be stopped after two gigabytes of memory 
had been exhausted. The first part of the table (Classical) 
suggests that the implementation does not behave well for 
large A^: the observed behaviour is exponential instead of 
polynomial. Even then, it is still much better than our im- 
plementation. Indeed, we have implemented only the case 
with rational values of m and for small N it often takes 
longer for our implementation to compute the result with 
this value than for the classical method to find the result 

1 All our tests have been run on an Intel Xeon at 3.6GHz. 
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Table 2: Timings (in sec.) for classical and compact 
versions of Zeilberger's algorithm 

with a formal m. However, things change as N gets larger: 
the predicted behaviour is well observed. When N is mul- 
tiplied by 2, the time is multiplied by slightly more than 2. 
Had we implemented the baby-step / giant-step version of bi- 
nary splitting, the timings in the columns for random m 
would have been much better, since the time should be mul- 
tiplied by slightly more than \/2 from one line to the next. 
Our experiments with symbolic m show that so far, our com- 
plexity result is more of a theoretical nature: although the 
degrees of the coefficients of the equations grow like O(N), 
the constant in front of the O term is about 18 in this exam- 
ple, and a massive cancellation takes place in the final linear 
solving. The result has degrees that also grow like O(N), 
but with a much smaller constant, so that a direct resolution 
in 0(N 4 ) is much faster in this range than our <D(N 2 ). 
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